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Abstract 

In  this  paper,  we  further  analyze,  test,  modify  and  improve  the  high  order  WENO 
(weighted  essentially  non-oscillatory)  finite  difference  schemes  of  Liu,  Osher  and  Chan 
[9].  It  was  shown  by  Liu  et  al.  that  WENO  schemes  constructed  from  the  r*'*  order 
(in  norm)  ENO  schemes  are  (r+l)^^  order  accurate.  We  propose  a  new  way  of 
measuring  the  smoothness  of  a  numerical  solution,  emulating  the  idea  of  minimizing 
the  total  variation  of  the  approximation,  which  results  in  a  5*^  order  WENO  scheme 
for  the  case  r  =  3,  instead  of  the  4^^  order  with  the  original  smoothness  measurement 
by  Liu  et  al.  This  5*^  order  WENO  scheme  is  as  fast  as  the  4*^  order  WENO  scheme 
of  Liu  et  al.  and,  both  schemes  are  about  twice  as  fast  as  the  4^^  order  ENO  schemes 
on  vector  supercomputers  and  as  fast  on  serial  and  parallel  computers.  For  Euler  sys¬ 
tems  of  gas  dynamics,  we  suggest  to  compute  the  weights  from  pressure  and  entropy 
instead  of  the  characteristic  values  to  simplify  the  costly  characteristic  procedure.  The 
resulting  WENO  schemes  are  about  twice  as  fast  as  the  WENO  schemes  using  the  char¬ 
acteristic  decompositions  to  compute  weights,  and  work  well  for  problems  which  do  not 
contain  strong  shocks  or  strong  reflected  waves.  We  also  prove  that,  for  conservation 
laws  with  smooth  solutions,  all  WENO  schemes  are  convergent.  Many  numerical  tests, 
including  the  ID  steady  state  nozzle  flow  problem  and  2D  shock  entropy  wave  inter¬ 
action  problem,  are  presented  to  demonstrate  the  remarkable  capability  of  the  WENO 
schemes,  especially  the  WENO  scheme  using  the  new  smoothness  measurement,  in  re¬ 
solving  complicated  shock  and  flow  structures.  We  have  also  applied  Yang’s  artificial 
compression  method  to  the  WENO  schemes  to  sharpen  contact  discontinuities. 
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1  Introduction 


In  this  paper,  we  further  analyze,  test,  modify  and  improve  the  WENO  (weighted  essentially 
non-oscillatory)  finite  difference  schemes  of  Liu,  Osher  and  Chan  [9]  for  the  approximation 
of  hyperbolic  conservation  laws  of  the  type: 


Ut  +  divf(u)  =  0 


(1.1) 


or  perhaps  with  a  forcing  term  g(u, x,  t)  on  the  right  hand  side.  Here  u  =  (wi, . . . ,  Um), f  = 
(fi,...,fd),x=  (a;i,...,x<i)  and  t  >  0. 

WENO  schemes  are  based  on  ENO  (essentially  non-oscillatory)  schemes,  which  were  first 
introduced  by  Harten,  Osher,  Engquist  and  Chakravarthy  [5]  in  the  form  of  cell  averages. 
The  key  idea  of  ENO  schemes  is  to  use  the  “smoothest”  stencil  among  several  candidates  to 
approximate  the  fluxes  at  cell  boundaries  to  a  high  order  accuracy  and  at  the  same  time  to 
avoid  spurious  oscillations  near  shocks.  The  ceU-averaged  version  of  ENO  schemes  involves 
a  procedure  of  reconstructing  point  values  from  cell  averages  and  could  become  complicated 
and  costly  for  multi-dimensional  problems.  Later,  Shu  and  Osher  [14,  15]  developed  the 
flux  version  of  ENO  schemes  which  does  not  require  such  a  reconstruction  procedure.  We 
win  formulate  the  WENO  schemes  based  on  this  flux  version  of  ENO  schemes.  The  WENO 
schemes  of  Liu  et  al.  [9]  are  based  on  the  ceU  averaged  version  of  ENO  schemes. 

For  applications  involving  shocks,  second  order  schemes  are  usually  adequate  if  only 
relatively  simple  structures  are  present  in  the  smooth  part  of  the  solution  (e.g.  the  shock 
tube  problem).  However,  if  a  problem  contains  rich  structures  as  well  as  shocks,  (e.g. 
the  shock  entropy  wave  interaction  problem  in  Example  4,  Section  8.3),  high  order  shock 
capturing  schemes  (order  of  at  least  three)  are  more  efficient  than  low  order  schemes  in 
terms  of  CPU  time  and  memory  reqidrements. 

ENO  schemes  are  uniformly  high  order  accurate  right  up  to  the  shock  and  are  very 
robust  to  use.  However,  they  also  have  certain  drawbacks.  One  problem  is  with  the  freely 
adaptive  stencil,  which  could  change  even  by  a  round-off  error  perturbation  near  zeroes  of 
the  solution  and  its  derivatives.  Also,  this  free  adaptation  of  stencils  is  not  necessary  in 
regions  where  the  solution  is  smooth.  Another  problem  is  that  ENO  schemes  are  not  cost 
effective  on  vector  supercomputers  such  as  the  CRAY  C-90  because  the  stencil  choosing 
step  involves  heavy  usage  of  logical  statements,  which  perform  poorly  on  such  machines. 
The  first  problem  could  reduce  the  accuracy  of  ENO  schemes  for  certain  functions  [12], 
however  this  can  be  remedied  by  embedding  certain  parameters  (e.g.  threshold  and  biasing 
factor)  into  the  stencil  choosing  step  so  that  the  preferred  linearly  stable  stencil  is  used  in 
regions  away  from  discontinuities.  See  [1,  3,  13]. 

WENO  scheme  of  Liu,  Osher  and  Chan  [9]  is  another  way  to  overcome  these  drawbacks 
while  keeping  the  robustness  and  high  order  accuracy  of  ENO  schemes.  The  idea  is  the 
following:  instead  of  approximating  the  numerical  flux  using  only  one  of  the  candidate 


stencils,  one  uses  a  convex  combination  of  all  the  candidate  stencils.  Each  of  the  candidate  ^  _ 

stencils  is  assigned  a  weight  which  determines  the  contribution  of  this  stencil  to  the  final 
approximation  of  the  numerical  flux.  The  weights  can  be  defined  in  such  a  way  that  in 
smooth  regions  it  approaches  certain  optimal  weights  to  achieve  a  higher  order  of  accuracy 
(a  order  ENO  scheme  leads  to  a  (2r-l)*^  order  WENO  scheme  in  the  optimal  case),  while 
in  regions  near  discontinuities,  the  stencils  which  contain  the  discontinuities  are  assigned  a 
nearly  zero  weight.  Thus,  the  essentially  non-oscillatory  property  is  achieved  by  emulating 

ENO  schemes  around  discontinrdties  and  a  higher  order  of  accuracy  is  obtained  by  emulating  — 

. . ...jilitiy  (^des 

lAvali  acd/or 

1  Mst,  j  Specie 


□  □ 


upstream  central  schemes  with  the  optimal  weights  away  from  discontinuities.  WENO 
schemes  completely  remove  the  logical  statements  that  appear  in  the  ENO  stencil  choosing 
step.  As  a  result,  the  WENO  schemes  run  at  leeist  twice  as  fast  as  ENO  schemes  (see 
Section  7)  on  vector  machines  (e.g.  CRAY  C-90)  and  are  not  sensitive  to  round-off  errors 
that  arise  in  actual  computation.  Atkins  [1]  also  has  a  version  of  ENO  schemes  nsing  a 
different  weighted  average  of  stencils. 

Another  advantage  of  WENO  schemes  is  that  its  flux  is  smoother  than  that  of  the  ENO 
schemes.  This  smoothness  enables  us  to  prove  convergence  of  WENO  schemes  for  smooth 
solutions  using  Strang’s  technique  [18],  see  Section  6.  According  to  our  numerical  tests,  this 
smoothness  also  helps  the  steady  state  calculations,  see  Example  4  in  Section  8.2. 

In  [9],  the  order  of  accuracy  shown  in  the  error  tables  (Table  1-5  in  [9])  seemed  to  suggest 
that  the  WENO  schemes  of  Liu  et  al.  are  more  accurate  than  what  the  truncation  error 
analysis  indicated.  In  Section  2,  we  carry  out  a  more  detailed  error  analysis  for  the  WENO 
schemes  and  find  that  this  “super-convergence”  is  indeed  superficial:  the  “higher”  order  is 
caused  by  larger  error  on  the  coarser  grids  instead  of  smaller  error  on  the  finer  grids.  Our 
error  analysis  also  suggests  that  the  WENO  schemes  can  be  made  more  accurate  than  those 
in  [9]. 

Since  the  weight  on  a  candidate  stencil  has  to  vary  according  to  the  relative  smoothness 
of  this  stencil  to  the  other  candidate  stencils,  the  way  of  evaluating  the  smoothness  of  a 
stencil  is  crucial  in  the  definition  of  the  weight.  In  Section  3,  we  introduce  a  new  way  of 
measuring  the  smoothness  of  the  numerical  solution  which  is  based  upon  mmimiying  the 
norm  of  the  derivatives  of  the  reconstruction  polynomials,  emulating  the  idea  of  TninimiTiiTig 
the  total  variation  of  the  approximations.  This  new  measurement  gives  the  optimal  5*^ 
order  accurate  WENO  scheme  when  r  =  3  (the  smoothness  measurement  in  [9]  gives  a  4*^ 
order  accurate  WENO  scheme  for  r  =  3). 

Although  the  WENO  schemes  are  faster  than  ENO  schemes  on  vector  supercomputers, 
they  are  only  as  fast  as  ENO  schemes  on  serial  computers.  In  Section  4,  we  present  a  simpler 
way  of  computing  the  weights  for  the  approximation  of  Euler  systems  of  gas  dynamics.  The 
simplification  is  aimed  at  reducing  the  floating  point  operations  in  the  costly  but  necessary 
characteristic  procedure  and  is  motivated  by  the  following  observation:  the  only  nonlinearity 
of  a  WENO  scheme  is  in  the  computation  of  the  weights.  We  suggest  the  use  of  pressure 
and  entropy  to  compute  the  weights  instead  of  the  local  characteristic  quantities.  In  this 
way  one  can  exploit  the  linearity  of  the  rest  of  the  scheme.  The  resulting  WENO  scheme 
(r  =  3)  is  about  twice  as  fast  as  the  original  WENO  scheme  which  uses  local  characteristic 
quantities  to  compute  the  weights  (see  Section  7).  The  same  idea  can  also  be  applied  to 
the  original  ENO  schemes.  Namely,  we  can  use  the  undivided  differences  of  pressure  and 
entropy  to  replace  the  local  characteristic  quantities  to  choose  the  ENO  stencil.  This  has 
been  tested  numerically  but  the  results  are  not  included  in  this  paper  since  the  main  topic 
here  is  the  WENO  schemes. 

WENO  schemes  have  the  same  smearing  at  contact  discontinuities  as  ENO  schemes. 
There  are  mainly  two  techniques  for  sharpening  the  contact  discontinuities  for  ENO  schemes. 
One  is  Harten’s  subcell  resolution  [4]  and  the  other  is  Yang’s  artificial  compression  (slope 
modification)  [20].  Both  were  introduced  in  the  cell  average  context.  Later,  Shu  and 
Osher  [15]  translated  them  into  the  point  value  framework.  In  one  dimensional  problems, 
subcell  resolution  technique  works  slightly  better  than  the  artificial  compression  method. 
However,  for  two  or  higher  dimensional  problems,  the  latter  is  found  to  be  more  effective  and 
easier  to  use  [15].  We  will  highlight  the  key  procedures  of  applying  the  artificial  compression 
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method  to  the  WENO  schemes  in  Section  5. 

In  Section  8,  we  test  the  WENO  schemes  (both  the  WENO  schemes  of  Liu  et  al.  and 
the  modified  WENO  schemes),  on  several  ID  and  2D  model  problems  and  compare  them 
with  ENO  schemes  to  examine  their  capability  in  resolving  shock  and  complicated  flow 
structures. 

We  conclude  this  paper  by  a  brief  summary  in  Section  9. 

The  time  discretization  of  WENO  schemes  will  be  implemented  by  a  class  of  high  order 
TVD  Runge-Kutta  type  methods  developed  by  Shu  and  Osher  [14].  To  solve  the  following 
ordinary  differential  equation: 

du  ^ 

Tt  =  (1-2) 

where  L{u)  is  a  discretization  of  the  spatial  operator,  the  third  order  TVD  Runge-Kutta  is 
simply: 


u(i)  =  tt” -t- AtL(u") 

+  l„(i)  +  lAtX(«W)  (1,3) 

+  |„(2)  +  BAti(«(2)) 

o  o  o 

Another  useful,  although  not  TVD,  fourth  order  Runge-Kutta  scheme  is: 
uW  =  u^  +  ^ML(u^) 

u(2)  = 

=  u”  -I-  AtL{u^‘^^)  (1.4) 

^  -1-  -t-  2u^^'>  +  -t-  ^AtL(u^^^) 

This  fourth  order  Rimge-Kutta  scheme  can  be  made  TVD  by  an  increa.se  of  operation 
counts  [14].  We  wfll  mainly  use  these  two  Runge-Kutta  schemes  in  our  numerical  tests  in 
Section  8.  The  third  order  TVD  Runge-Kutta  scheme  wfll  be  referred  to  as  “RK-S”  while 
the  fourth  order  (non-TVD)  Runge-Kutta  scheme  will  be  referred  to  as  “RK-4”. 

2  The  WENO  Schemes  of  Liu,  Osher  and  Chan 

In  this  section,  we  use  the  flux  version  of  ENO  schemes  as  our  basis  to  formulate  WENO 
schemes  of  Liu  et  al.  and  analyze  their  accuracy  in  a  different  way  from  that  used  in  [9]. 
We  use  one  dimensional  scalar  conservation  laws  (i.e.  d=m=l  in  (1-1)  )  as  an  example: 

+  /(«)ar  =  0  (2.1) 

Let  us  discretize  the  space  into  uniform  intervals  of  size  Ax  and  denote  Xj  =  jAx. 
Various  quantities  at  xj  will  be  identified  by  the  subscript  j.  The  spatial  operator  of  the 
WENO  schemes,  which  approximates  —f{u)x  at  Xj,  will  take  the  following  conservative 
form: 

^  ^  “  Az  ~ 
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where  the  numerical  fliix  /i+l  approximates  =  h{Xj^x)  to  a  high  order  with  h{x) 
implicitly  defined  by  [15] 

1  rx+Aa:/2 

(2.3) 

We  can  actually  assume  f'{u)  >  0  for  all  u  in  the  range  of  our  interest.  For  a  general 
flux,  i.e.  f{y)  0,  one  can  split  it  into  two  parts  either  globally  or  locally: 

f(u)  =  f+{u)  +  f-  (u)  (2.4) 

where  >  0  and  — ^  <  0.  For  example,  one  can  define 

/*(«)  =  \  (/(«)  ±  (2.5) 

where  a  =  max|/'(u)|  and  the  maximum  is  taken  over  the  whole  relevant  range  of  u.  This 
is  the  global  Lax- Friedrichs  (LF)  flux  splitting.  For  other  flux  splittings,  especially  the  Roe 
flux  splitting  with  entropy  fix  (RF),  see  [15]  for  details.  Let  /+ i  and  /“  i  be,  resp.  the 

numerical  fluxes  obtained  from  the  positive  and  negative  parts  of  f{u),  we  then  have: 

/i+i  =  //tj  +  4.J  (2-8) 

Here  we  will  only  describe  how  /^i  is  computed  in  [9]  on  the  basis  of  the  flux  version 

of  ENO  schemes.  For  simplicity,  we  will  drop  the  “-I-”  sign  in  the  superscript.  The  formulas 
for  the  negative  part  of  the  split  flux  are  symmetric  (with  respect  to  Xj+i )  and  will  not  be 
shown. 

As  we  know,  the  order  (in  sense)  ENO  scheme  chooses  one  “smoothest”  stencil 
from  r  candidate  stencils  and  uses  only  the  chosen  stencil  to  approximate  the  flux  h-,i. 

3  r  2 

Let’s  denote  the  r  candidate  stencils  by  5*,  fc  =  0, 1, . . . ,  r- 1  where 

Sk  —  (®j-|-fc-r+l )  r+25  ’ '  ‘  i  ®j+fc) 

If  the  Stencil  5*  happens  to  be  chosen  as  the  ENO  interpolation  stencil,  then  the  order 
ENO  approximation  of  h-j^i  is: 


where 


/j+i  —  qk{fj+k-r+l ,  •  -  • ,  fj+k) 


9ki9o,  •  •  • ,  9t-i)  =  ^  al  igi 


Here  /,0  <k,l  <  r-1  are  constant  coefl5cients.  For  later  use,  we  provide  these  coefficients 
for  r  =  2,3  in  Table  1. 

To  just  use  the  one  smoothest  stencil  among  the  r  candidates  for  the  approximation 
of  ,  is  very  desirable  near  discontinuities  because  it  prohibits  the  usage  of  information 
on  discontinuous  stencUs.  However,  it  is  not  so  desirable  in  smooth  regions  because  all  the 
candidate  stencils  carry  equally  smooth  information  and  thus  can  be  used  together  to  give 
a  higher  order  (higher  than  r,  the  order  of  the  base  ENO  scheme)  approximation  to  the  flux 
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Table  1:  Coefficients 


r 

k 

/  =  0 

1  =  1 

1  =  2 

2 

0 

-1/2 

3/2 

1 

1/2 

1/2 

3 

0 

1/3 

-7/6 

11/6 

1 

-1/6 

5/6 

1/3 

2 

1/3 

5/6 

-1/6 

In  fact,  one  could  use  all  the  r  candidate  stencils,  which  all  together  contain  (2r— 1) 
grid  values  of  /  to  give  a  (2r— 1)*^  order  approximation  of 

/j+i  =  9r-l^(/i-r+l  >  •  •  • ,  /i+r-l)  (2.9) 

which  is  just  the  numerical  flux  of  a  (2r--l)*'‘  order  upstream  central  scheme.  As  we  know, 
high  order  upstream  central  schemes  (in  space)  combined  with  high  order  Runge-Kutta 
methods  (in  time),  are  stable  and  dissipative  under  appropriate  CFL  numbers  and  thus 
are  convergent,  according  to  Strang’s  convergence  theory  [18]  when  the  solution  of  (1.1)  is 
smooth  (see  Section  6).  The  above  facts  suggest  that,  one  could  use  the  (2r  — 1)*^  order 
upstream  central  scheme  in  smooth  regions  and  only  use  the  order  ENO  scheme  near 
discontinuities. 

As  in  (2.7),  each  of  the  stencils  can  render  an  approximation  of  h-,i.  If  the  stencil  is 

smooth,  this  approximation  is  order  accurate,  otherwise,  it  is  less  accurate  or  even  not 
accurate  at  all  if  the  stencil  contains  a  discontinuity.  One  could  assign  a  weight  Uk  to  each 
candidate  stencil  Sk^  k  =  0, l,.,,,r  — 1  and  use  these  weights  to  combine  the  r  different 
approximations  to  obtain  the  final  approximation  of  /i- .  i  as  follows: 

r—l 

fj+l  ~  ‘^kiUfj+k-r+l ,  •  •  • ,  fj+k)  (2.10) 

k=0 

where  ql{fj+k-r+i,”-,fj+k)  is  defined  in  (2.8).  To  achieve  the  essentially  non-osciUatory 
property,  one  then  requires  the  weights  to  adapt  to  the  relative  smoothness  of  /  on  each 
candidate  stencil  such  that  any  discontinuous  stencil  is  effectively  assigned  a  zero  weight.  In 
smooth  regions,  one  can  adjust  the  weight  distribution  such  that  the  resxdting  approximation 
of  the  flux  fj_^i  is  as  close  as  possible  to  that  given  in  (2.9). 

Simple  algebra  gives  the  coefficients  Cl  such  that 

r—l 

9r-l  ifj-r+1  ? '  ‘  ?  /j+r-l)  ~  ?  *  *  * ;  f.i+k)  (2.11) 

fc=0 

and  ^k  —  ^  r  >  2.  For  r  =  2, 3,  these  coefficients  are  given  in  Table  2. 

Comparing  (2.11)  with  (2.10),  we  get: 

r  —  l 

fj+k  =  Qr-l^(fj-r+l'>‘  ”  1  fj+T-l)  +  ^(‘^i  —Cl)ql{fj.kk-r+l,---'>fj+k)  (2.12) 

fc=0 
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Table  2:  Optimal  weights  C^. 


Cl 

k=0 

k=l 

k=2 

r=2 

1/3 

2/3 

— 

r=3 

1/10 

6/10 

3/10 

Recalling  (2.9),  we  see  that,  the  first  term  on  the  right  hand  side  of  the  above  equation  is  a 
(2r-l)‘^  order  approximation  of  Since  C"!  =  1,  if  we  require  E*=o^Jfc  =  1,  the 
last  summation  term  can  be  written  as 


-  Cl)(gl(fj+k-r+l ,  •  •  • ,  fj+k)  - 

k=0  * 

Each  term  in  the  last  summation  can  be  made  0(h^^~^)  if 

wjt  =  q  +  0(h^-^) 


(2.13) 


(2.14) 


OJk  =  - 

ao  H - 1-  Q!r-i 

(2.15) 

1! 

o 

1 

(2.16) 

for  A:  —  0, 1,  —  Here,  k  =  Ax.  Thus  will  bear  the  name  of  optimal  weight. 

The  question  now  is  how  to  define  the  weight  such  that  (2.14)  is  satisfied  in  smooth 
regions  while  essentially  non-oscillatory  property  is  achieved.  In  [9],  the  weight  for 
stencil  Sk  is  defined  by 


where 


Here  e  is  a  positive  real  number  which  is  introduced  to  avoid  the  denominator  to  become 
zero  ( in  our  later  tests,  we  will  take  e  =  10“®  );  the  power  p  will  be  discussed  in  a  moment; 
ISk  in  (2.16)  is  a  smoothness  measurement  of  the  flux  function  on  the  candidate  stencil. 

It  is  easy  to  see  that  EI=o  =  1-  To  satisfy  (2.14),  it  sufiices  to  have  (through  a  Taylor 
expansion  analysis): 

ISk  =  D{1  +  O(fi’-i))  (2.17) 

for  A:  =  0, 1,  —  ,?*— 1  where  D  is  some  nonzero  quantity  independent  of  k. 

As  we  know,  an  ENO  scheme  chooses  the  “smoothest”  ENO  stencil  by  comparing  a 
hierarchy  of  undivided  differences.  This  is  because  these  undivided  differences  can  be  used 
to  measure  the  smoothness  of  the  numerical  flux  on  a  stencil.  In  [9],  ISk  is  defined  as 


1=1  .=1  ’*  “  ^ 

where  /[*,  •]  is  the  undivided  difference: 

/[i,0]  =  fj 

f[jJ]  =  k  = 

For  example,  when  r  =  2,  we  have 

ISk  =  if\j  +  k-l,l])^  k  =  0, 


(2.18) 


A:  =  l,...,r-l. 


(2.19) 
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When  r  =  3,  (2.18)  gives 

ISk  =  I  {if[j  +  k-2,  l]f  +  ifU  +  k-l,  1])2)  +  {f[j  +  k-2, 2]f  k  =  0,1,2  (2.20) 

In  smooth  regions,  Taylor  expansion  analysis  of  (2.18)  gives 

ISk  =  {fh)\l  +  0(h))  k  =  0,l,...,r-l.  (2.21) 

where  f  =  f\uj).  Note  the  0{h)  term  is  not  0{h''~^)  that  we  would  want  to  have  (see 
(2.17)).  Thus  in  smooth  monotone  regions,  i.e.  f  ^  0,  we  have: 

Uk  =  Cl-\-0{h)  fc  =  0,l,...,r-l.  (2.22) 

Recalling  (2.12),  we  see  that  the  WENO  schemes  with  the  smoothness  mea.surement  given 
by  (2.18)  is  (r  +  1)*^  order  accurate  in  smooth  monotone  regions  of  f{u{x)).  This  result 
was  proven  in  [9]  using  a  different  approach.  For  r  =  2,  this  is  optimal  in  the  sense  that 
the  3’^'^  order  upstream  central  scheme  is  approximated  in  most  smooth  regions.  However, 
this  is  not  optimal  for  r  =  3,  for  which  this  measurement  can  only  give  4*^  order  accuracy 
while  the  optimal  upstream  central  scheme  is  5*^  order  accurate.  We  will  introduce  a  new 
measurement  in  the  next  section  which  will  residt  in  an  optimal  order  accurate  WENO 
scheme  for  the  r  =  3  case. 

When  r  =  3,  Taylor  expansion  of  (2.20)  gives: 

ISo  =  ^((A-^/'V)2  +  (/h-l/'V)2)+(/"h2)2  +  C)(h5)  (2.23) 

ISx  =  \  ((A  -  \f"hy  +  {fh  +  +  {f'hy  +  0{h^)  (2.24) 

IS2  =  \  (ifh  +  i/'V)2  +  (fh  +  |/'V)2)  +  (fhy  +  0{h^)  (2.25) 

We  can  see  that  the  second  order  terms  are  different  from  stencil  to  stencil.  Thus  (2.22)  is 
no  longer  valid  at  critical  points  of  f(u(x))  which  implies  that  the  WENO  scheme  of  Liu 
et  al.  for  r  =  3  is  only  3’"'^  order  accurate  at  these  points.  In  fact,  the  weights  computed 
from  the  smoothness  measurement  (2.18)  diverge  far  away  from  the  optimal  weights  near 
critical  points  (see  Figure  1  in  the  next  section)  on  coarse  grids  (10  to  80  grid  points  per 
wave).  But  on  fine  grids,  since  the  smoothness  measurements  ISk  for  all  k  are  relatively 
smaller  than  the  non-zero  constant  e  in  (2.16),  the  weights  become  close  to  the  optimal 
weights.  Therefore  the  “super-convergence”  phenomena  which  appeared  in  Table  1-5  in  [9] 
are  caused  by  large  error  commitment  on  coarse  grids  and  less  error  commitment  on  finer 
grids  when  using  the  errors  of  the  5*^  order  central  scheme  as  reference  (see  Table  3  and  4). 

At  discontinuities,  it  is  typical  that  one  or  more  of  the  r  candidate  stencils  reside  in 
smooth  regions  of  the  numerical  solution  while  other  stencils  contain  the  discontinuities. 
The  size  of  the  discontinuities  is  always  0(1)  and  does  not  change  when  the  grid  is  refined. 
So  we  have  for  a  smooth  stencil  Sk, 

ISk  =  0(h^P)  (2.26) 

and  for  a  non-smooth  stencil  Si, 

ISi  =  0(1)  (2.27) 


7 


From  the  definition  of  the  weights  (2.15),  we  can  see  that,  for  this  non-smooth  stencil  Si, 
the  corresponding  weight  w/  satisfies 


wz  =  0(/i2P)  (2.28) 

Therefore  for  small  h  and  any  positive  integer  power  p,  the  weight  assigned  to  the  non¬ 
smooth  stencil  vanishes  as  h  — >  0.  Note,  if  there  is  more  than  one  smooth  stencil  in  the 
r  candidates,  from  the  definition  of  the  weights  in  (2.15),  we  expect  each  of  the  smooth 
stencils  will  get  a  weight  which  is  0(1).  In  this  case,  the  weights  do  not  exactly  resemble 
the  “ENO  digital  weights”.  However,  if  a  stencil  is  smooth,  the  information  that  it  contains 
is  useful  and  should  be  utilized.  In  fact,  in  our  extensive  numerical  experiments,  we  find 
the  WENO  schemes  in  [9]  work  very  well  at  shocks.  We  also  find  that  p  =  2  is  adequate 
to  obtain  essentially  non-osciUatory  approximations  at  least  for  r  =  2,3,  although  it  is 
suggested  in  [9]  that  p  should  be  taken  as  r,  the  order  of  the  base  ENO  schemes.  We  will 
use  p  =  2  for  all  our  numerical  tests. 

In  summary,  WENO  schemes  of  Liu  et  al.  defined  by  (2.10),  (2.15)  and  (2.18)  have  the 
following  properties: 

1.  They  involve  no  logical  statements  which  appear  in  the  base  ENO  schemes. 

2.  The  WENO  scheme  based  on  the  order  ENO  scheme  is  (r-|-l)*^  order  accurate 
in  smooth  monotone  regions,  although  this  is  still  not  as  good  as  the  optimal  order 
(2r-l)‘'‘. 

3.  They  achieve  the  essentially  non-osciUatory  property  by  emulating  ENO  schemes  at 
discontinuities. 

4.  They  are  smooth  in  the  sense  that  the  numerical  fliLX  /• ,  i  is  a  smooth  function  of  all 

J't  2 

its  arguments  (For  a  general  flux,  this  is  also  true  if  a  smooth  flux  splitting  method 
is  used,  e.g.  global  Lax-Friedrichs  flux  splitting). 


3  A  New  Smoothness  Measurement 

In  this  section,  we  present  a  new  way  of  measuring  the  smoothness  of  the  numerical  solution 
on  a  stencil  which  can  be  used  to  replace  (2.18)  to  form  a  new  weight. 

As  we  know,  on  each  stencil  Sk,  we  can  construct  a  (r— 1)*^  order  interpolation  poly¬ 
nomial,  which  if  evaluated  at  x  =  renders  the  approximation  of  given  in  (2.7). 

Since  total  variation  is  a  good  measurement  for  smoothness,  it  would  be  desirable  to  mini¬ 
mize  the  total  variation  for  the  approximation.  Consideration  of  a  smooth  flux  and  of  the 
role  of  higher  order  variations  leads  us  to  the  following  measurement  for  smoothness:  let 
the  interpolation  polynomial  on  stencil  Sk  be  qk(x),  we  define 

(3.1) 

1=1 

where  is  the  derivative  of  qk(x).  The  right  hand  side  of  (3.1)  is  just  a  sum  of 
the  norms  of  all  the  derivatives  of  the  interpolation  polynomial  qk{x')  over  the  interval 
The  term  0~'‘-  is  to  remove  h-dependent  factors  in  the  derivatives  of  the 
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polynomials.  This  is  similar  to  but  smoother  than  the  total  variation  measurement  based 
on  the  norm.  It  also  renders  a  more  accurate  WENO  scheme  for  the  case  r  =  3,  when 
used  with  (2.15)  and  (2.16). 

When  r  =  2,  (3.1)  gives  the  same  measurement  as  (2.18).  However,  they  become 


different  for  r  >  3.  For  r  =  3,  (3.1)  gives 

150  -  + +  (3.2) 

151  =  (3.3) 

ISz  =  j|(/i-2/,+.+/j+2)^  +  i(3/,-4/,+,  +  /,+,)=  (3.4) 

In  smooth  regions.  Taylor  expansion  of  (3.2)-(3.4)  gives,  resp. 

150  =  +  \(2f'h  -  +  0(h^)  (3.5) 

151  =  +  \{2fh  +  +  0{h^)  (3.6) 

152  =  ^U"hy  +  \{2rh  -  \f"hy  +  0{h^)  (3.7) 

where  f"  =  If  ^  0,  then 

ISk  =  (fh)\l  +  0(/i2))  fc  =  0, 1,2  (3.8) 


which  means  the  weights  (see  (2.15))  resulting  from  this  measurement  satisfy  (2.17)  for 
r  =  3,  thus  we  obtain  a  5^^  order  (the  optimal  order  for  r  =  3)  accurate  WENO  scheme. 

Moreover,  this  measurement  is  also  more  accurate  at  critical  points  of  f{u{x)).  When 
/'  =  0,  we  have 

ISk  =  ||(/'V)2(1  +  0{h^))  k  =  0,1,2  (3.9) 

which  implies  that  the  weights  resulting  from  the  measurement  (3.1)  are  also  5*^  order 
accurate  at  critical  points. 

To  illustrate  the  different  behavior  of  the  two  measurements  (i.e.  (2.18)  and  (3.1))  for 
r  =  3  in  smooth  monotone  regions,  near  critical  points  or  near  discontinuities,  we  compute 
the  weights  uq,ui  and  a;2  for  the  following  function: 

,  f  sin27ra:j  if  0  <  xj  <  0.5,  .  , 

\  1  -  sin27rxj  if  0.5  <  Xj  <1.  \  •  ) 

at  all  half  grid  points  where  Xj  =  jAx,  =  Xj  +  Ax/2  and  Ax  =  1/40.  We 
display  the  weights  a;o  and  ui  in  Figure  1.  (a;2  =  1  —  wq  —  a;i  is  omitted  in  the  picture). 
Note  the  optimal  weight  for  uq  is  Cq  =  0.1  and  for  wi  is  Cf  =  0.6.  We  can  see  that  the 
weights  computed  with  (2.18)  (referred  to  as  the  original  measurement  in  Figure  1)  are 
far  less  optimal  than  those  with  the  new  measurement  especially  around  the  critical  points 
X  =  |,|.  However,  near  the  discontinuity  x  =  ^,  the  two  measurements  behave  similarly: 
the  discontinuous  stencil  always  gets  an  almost  zero  weight.  Moreover,  for  the  grid  point 
immediately  to  the  left  of  the  discontinuity,  ljq  ^  j  and  which  means,  when  only  one 

of  the  three  stencils  is  non-smooth,  the  other  two  stencils  get  0(1)  weights.  Unfortunately, 
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® - Original  Measurement  (o^ 

0  Original  Measurement  o, 

. .  New  Measurement 


Figure  1:  A  comparison  of  the  two  smoothness  measurements. 

these  weights  do  not  approximate  a  4*^  order  scheme  at  this  point.  A  similar  situation 
happens  to  the  point  just  to  the  right  of  the  discontinuity. 

For  simplicity  of  notations,  we  use  WENO-X-3  to  stand  for  the  3’"'^  order  WENO  scheme 
(i.e.  r  =  2,  for  which  the  original  and  new  smoothness  measurement  coincide)  where  X=LF, 
Roe,  RF  refers  resp.  to  the  global  Lax-Friedrichs  flux  splitting.  Roe’s  flux  splitting  and  Roe’s 
flux  splitting  with  entropy  fix;  The  accuracy  of  this  scheme  has  been  tested  in  [9].  We  will 
use  WENO-X-4  to  represent  the  4**'  order  WENO  scheme  of  Liu  et  al.  (i.e.  r  =  3  with 
the  original  smoothness  measurement  of  Liu  et  al.)  and  WENO-X-5  to  stand  for  the  5*^ 
order  WENO  scheme  restilting  from  the  new  smoothness  measurement.  In  later  sections, 
we  will  also  use  ENO-X-Y  to  denote  conventional  ENO  schemes  of  “Y”th  order  with  “X” 
flux  splitting.  We  caution  the  reader  that  the  orders  here  are  in  sense.  So  ENO-RF-4  in 
our  notation  refers  to  the  same  scheme  as  ENO-RF-3  in  [15]. 

In  the  following  we  test  the  accuracy  of  WENO  schemes  on  the  linear  equation: 

Uf  +  Ux  =  0  -l<a:<l  (3.11) 

u(a:,0)  =  uo{x)  periodic.  (3.12) 

In  Table  3,  we  show  the  errors  of  the  two  schemes  at  i  =  1  for  the  initial  condition  uo(x)  = 
sin(xa;)  and  compare  them  with  the  errors  of  the  5*^  order  upstream  central  scheme  (referred 
to  as  CENTRAL-5  in  the  following  tables).  We  can  see  that  WENO-RF-4  is  more  accurate 
than  WENO-RF-5  on  the  coarsest  grid  (N=10)  but  becomes  less  accurate  than  WENO-RF- 
5  on  the  finer  grids.  Moreover,  WENO-RF-5  gives  the  expected  order  of  accuracy  starting 
at  about  40  grid  points.  In  this  example  and  the  one  for  Table  4,  we  have  adjusted  the  time 
step  to  At  ~  (Aa;)!  so  that  the  4*^  order  Runge-Kutta  in  time  is  effectively  5*^  order. 

In  Table  4,  we  show  errors  for  the  initial  condition  uo(®)  =  sin^(?r®).  Again  we  see  that 
WENO-RF-4  is  more  accurate  than  WENO-RF-5  on  the  coarsest  grid  (N=20)  but  becomes 
less  accurate  than  WENO-RF-5  on  finer  grids.  The  order  of  accuracy  for  WENO  settles 
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Table  3:  Accuracy  on  ut  +  Uj.  =  0  with  uo(®)  =  sm(7rx). 


Method 

N 

Loo  error 

Loo  order 

Li  error 

Li  order 

■El 

1.31e-2 

- 

7.93e-3 

- 

El 

3.00e-3 

2.13 

1.32e-3 

2.59 

WENO-RF-4 

Bl 

4.27e-4 

2.81 

1.56e-4 

3.08 

80 

5.17e-5 

3.05 

1.13e-5 

3.79 

160 

4.99e-6 

3.37 

6.88e-7 

4.04 

320 

3.44e-7 

3.86 

2.74e-8 

4.65 

iSI 

2.98e-2 

- 

1.60e-2 

- 

El 

1.45e-3 

4.36 

7.41e-4 

4.43 

WENO-RF-5 

m 

4.58e-5 

4.99 

2.22e-5 

5.06 

80 

1.48e-6 

4.95 

6.91e-7 

5.01 

160 

4.41e-8 

5.07 

2.17e-8 

4.99 

320 

1.35e-9 

5.03 

6.79e-10 

5.00 

■El 

4.98e-3 

- 

3.07e-3 

- 

1.60e-4 

4.96 

9.92e-5 

4.95 

CENTRAL-5 

n 

5.03e-6 

4.99 

3.14e-6 

4.98 

80 

1.57e-7 

5.00 

9.90e-8 

4.99 

160 

4.91e-9 

5.00 

3.11e-9 

4.99 

320 

1.53e-10 

5.00 

9.73e-ll 

5.00 

down  later  than  in  the  previous  example.  Notice  that  this  is  the  example  for  which  ENO 
schemes  lose  their  accuracy  [12]. 

4  A  Simple  Way  for  Computing  Weights  for  Euler  Systems 

For  system  (1.1)  with  d  >  1,  the  derivatives  =  1,. .  .,d  are  approximated  dimension 

by  dimension:  for  example,  when  approximating  one  fixes  xt,l  >  1  and  uses  an  one 
dimensional  approximation  in  the  direction  of  a;i.  In  the  foEowing,  we  only  discuss  how  to 
approximate  ^  and  will  drop  the  index  “1”  for  simpEcity.  We  wiE  also  cissume  that  aE 

the  eigenvalues  of  the  Jacobian  ^  are  nonnegative  (a  condition  identical  to  /'  >  0  in  the 
scalar  equation).  For  a  general  flux,  one  can  spEt  it  locaEy  into  positive  and  negative  parts 
just  as  in  the  scalar  case.  The  formulas  for  the  negative  part  of  the  flux  WEI  be  omitted  due 
to  symmetry. 

For  systems  of  equations,  the  fluxes  are  usuaEy  approximated  in  the  (local)  char¬ 
acteristic  fields.  Let’s  take  to  be  some  average  Jacobian  at  e.g.,  the  arithmetic 

mean  * 

.  _df\ 

“  au  r=K+"i+i)/2  (4.1) 

or  for  Euler  systems,  the  Roe’s  mean  matrix  [11].  We  denote  by  r®  (colunm  vector)  and 
Is  (row  vector)  the  right  and  left  eigenvector  of  ,  resp.  Then  the  scalar  WENO 
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Table  4:  Accuracy  on  ttj  +  Uj;  =  0  with  uo{x)  =  sm‘*(5ra:). 


Method 

N 

Loo  error 

Loo  order 

Li  error 

Li  order 

20 

7.31e-2 

« 

3.29e-2 

- 

40 

2.48e-2 

1.56 

9.99e-3 

1.72 

WENO-RF-4 

80 

4.60e-3 

2.43 

1.44e-3 

2.79 

160 

3.59e-4 

3.68 

8.31e-5 

4.12 

320 

2.12e-5 

4.08 

3.06e-6 

4.76 

640 

1.51e-6 

3.81 

9.57e-8 

5.00 

20 

1.08e-l 

- 

4.91e-2 

- 

40 

8.90e-3 

3.60 

3.64e-3 

3.75 

WENO-RF-5 

80 

1.80e-3 

2.31 

5.00e-4 

2.86 

160 

1.22e-4 

3.88 

2.17e-5 

4.53 

320 

4.37e-6 

4.80 

6.17e-7 

5.14 

640 

9.79e-8 

5.48 

1.57e-8 

5.30 

20 

5.23e-2 

- 

3.35e-2 

- 

40 

2.47e-3 

4.40 

1.52e-3 

4.46 

CENTRAL-5 

80 

8.32e-5 

4.89 

5.09e-5 

4.90 

160 

2.65e-6 

4.97 

1.60e-6 

4.99 

320 

8.31e-8 

5.00 

4.99e-8 

5.00 

640 

2.60e-9 

5.00 

1.56e-9 

5.00 

scheme  can  be  applied  to  each  of  the  characteristic  fields.  For  example,  (2.10)  becomes 

r— 1 

^  ^  Is’fj+Ar)  (4*2) 

A:=0 

which  gives  the  numerical  flux  in  the  characteristic  field.  Here  0Jk,s^k  =  0,1,..., r-1 
are  the  weights  in  the  characteristic  field: 

r+l?  •  •  l)  (4'3) 

which  is  a  nonlinear  function,  (uk  is  defined  by  (2.15)).  The  numerical  fluxes  obtained  in 
each  characteristic  field  can  then  be  projected  back  to  the  component  space  by: 

m 

5=1 

Because  of  the  nonlinearity  of  the  weights  (see  (4.3)),  the  above  procedure  involves  many 
local  projections  (or  vector  vector  products).  In  fact,  these  projections  are  responsible 
for  most  of  the  floating  point  operations  of  WENO  schemes  (true  also  for  ENO  schemes). 
Moreover,  these  projections  can  not  be  avoided  if  the  weights  are  to  be  computed  from  the 
projected  quantities.  However,  if  the  weights  can  be  computed  from  other  quantities,  we 
then  can  exploit  the  linearity  of  the  rest  of  the  scheme  (e.g.  the  linearity  of  ql)  to  reduce 
the  number  of  floating  point  operations  because  the  only  nonlinear  part  of  WENO  schemes 
is  in  the  calculation  of  the  weights. 
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The  question  then  is  what  quantities  can  serve  as  replacements  of  the  projected  val¬ 
ues.  Obviously  for  each  characteristic  field,  the  replacing  quantity  must  indicate  the  jump 
discontinuities  in  that  field.  Although  such  quantities  are  yet  to  be  discovered  for  general 
systems  of  equations,  we  find,  after  an  extensive  searching  and  trial,  that  pressure  and  en¬ 
tropy  are  good  replacements  for  the  projected  values  when  Euler  systems  are  concerned,  at 
least  for  problems  without  strong  shocks  and  reflective  waves. 

Namely,  we  will  use  pressure  to  compute  the  weights  in  the  genuinely  nonlinear  charac¬ 
teristic  fields  (s  =  l,m)  and  use  entropy  for  the  linearly  degenerate  field(s)  (1  <  s  <  m). 
The  motivation:  (1).  The  pressure  does  not  jump  at  contact  discontinuities  but  always 
jumps  at  shocks;  (2).  The  entropy  jumps  at  contact  discontinuities  but  jumps  only  slightly 
at  a  weak  shock. 

Since  the  pressure  and  entropy  can  be  obtained  independent  of  the  characteristic  projec¬ 
tion  procedure,  we  can  reformulate  the  WENO  schemes  to  take  advantage  of  the  linearity 
of  the  rest  of  the  scheme.  Let’s  define 

r-l 

~  ^  ^  ^k,sQk(.^j+k—T+l  j  ■  ■  ■  >  ^j+k)  ^  —  1, . . . ,  m 

jt=0 

For  Euler  systems,  the  (1  <  s  <  m)  characteristic  field  is  linearly  degenerate.  These 
fields  have  the  same  characteristic  speed  (eigenvalue)  and  the  weights  are  all  computed  from 
the  entropy.  So  we  have  for  all  1  <  s  <  m: 

(^k,s=(^k,2  VA:  =  0,...,r-l. 

and  therefore 

=  ^+1.2 

for  all  1  <  s  <  m.  Combine  (4.2)  and  (4.4)  and  use  the  linearity  of  ql  to  take  out  l^,  we  get 

m 

5=1 

-  •^i+1,2))  '*1  +  (im  -  ^i+§,2))  Tm  +  -?V+1,2  (4-6) 

As  we  can  see,  we  only  need  two  projections  from  component  space  to  characteristic  space 
and  two  inverse  projections,  plus  the  few  operations  for  computing  55^  =  l,2,m. 

We  wiU  denote,  by  WEN0-LF-5-PS,  the  WENO  scheme  for  the  case  r  =  3,  which 
uses  pressure  and  entropy  for  weight  computation  in  conjunction  with  the  new  smoothness 
measurement  (3.1),  the  weights  (2.15)  and  global  Lax-Friedrichs  flux  splitting  (according 
to  our  numerical  tests,  the  original  smoothness  measurement  of  Liu  et  al.  does  not  perform 
well  at  shocks  when  combined  with  the  above  way  of  computing  weights). 

Accuracy  of  WENO-LF-4,  WENO-LF-5  and  WEN0-LF-5-PS  on  the  ID  Euler  system 
is  tested  using  an  initial  condition  which  produces  a  smooth  solution,  the  same  example 
used  in  Section  6.  The  result  is  similar  to  the  scalar  case  in  Table  3  and  thus  will  not  be 
shown. 
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Sh^u:pening  of  Contact  Discontinuities 


For  alinear,  constant  coefficient  problem  (  f{u)  =  au  in  (2.1) ),  Yang’s  artificial  compression 
method,  when  applied  to  the  WENO  schemes  is  simply  (  assuming  a  >  0  ): 

//+.  =  /rtl  +  S+l  (5.1) 

where  is  the  flux  obtained  by  one  of  the  methods  introduced  in  the  previous  three 
sections,  and 


s+l  =  ^ 


fj-i 


Here  /^i  is  obtained  by  the  same  method  for  /• ,  i  pretending  a  <  0;  m  is  the  usual 
minmod  function  defined  by 


^(®l  ?  *  ’  *  5  ®n)  — 


min  I  a,- 1,  if  s  =  $ign{a\)  =  •  •  •  =  $ign{an) 
otherwise; 


and  aj  is  given  by 

where  o:  is  a  positive  parameter.  We  will  use  a  =  33  as  suggested  by  Yang  [20]  in  all  our 
tests  in  Section  8,  though  this  parameter  can  be  tuned  to  optimize  the  results  for  individual 
problems.  The  case  of  c  <  0  can  be  treated  symmetrically  and  the  generalization  to  variable 
coefficient  or  nonlinear  problems  is  rather  straight  forward.  See  [15]  for  details. 

We  will  apply  the  above  sharpening  technique  only  to  contact  discontinuities  or  contact 
characteristic  field(s)  in  case  of  Euler  systems.  A  scheme  which  uses  the  above  artificial 
technique  will  be  denoted  by  adding  to  its  name  the  suflSx  “-A”,  e.g.  WENO-LF-5-A. 


6  Convergence  for  Smooth  Solutions 

As  we  can  see  from  the  previous  sections,  the  WENO  schemes  are  smooth  in  the  sense  that 
the  spatial  operator  L 

^  —  L{fj-r,fj-r+l,-  •  •,/j+r-l)  (6.1) 

is  infinitely  differentiable  to  any  of  its  arguments  (  see  (2.2),  (2.10),  (2.15),  (2.16)  and  (2.18) 
or  (3.1) ).  Here  r  >  2  is  the  Li  order  of  the  base  ENO  scheme.  In  case  of  a  general  flux,  if  a 
smooth  flux  splitting  is  used  (e.g.  the  global  Lax-Friedrichs  flux  splitting),  the  smoothness 
of  the  WENO  schemes  is  unchanged. 

Strang’s  theorem  (Theorem  I  in  [18])  implies  that,  for  a  conservation  law  whose  flux 
function  and  solution  have  enough  continuous  derivatives,  a  smooth,  consistent  scheme  is 
convergent  if  its  first  variation  (see  [18]  for  the  definition)  is  /2-stable. 

It  is  easy  to  see  that,  for  the  scalar  one  dimensional  conservation  law  (2.1)  with  f  >  0, 
the  spatial  operator  of  WENO  schemes  has  the  following  simple  first  variation  L 
i+’-i  Qi 

Izzj^r  * 

=  ‘  •  5  «i+r-l)  “  *  *  • ,  Uj+^-a))  (6.2) 
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because  •,%)  =  0  and  =  Cl  for  all  A:  =  0,...,r-  1  and  I  = 

j  -  +  r  -  1.  (6.2)  can  be  rewritten  into  a  summation  of  a  (2r  -  2)*^  order  central 

difference  and  a  (2r  —  1)*^  order  upwind  biased  difference. 

i  +  (6.3) 

where  ^  ®  is  the  (2r  —  1)*^  order  forward  difference  operator. 

Applying  the  classical  Fourier  analysis  to  the  first  variation,  we  see  that  the  (2r  —  2)*^  order 
central  difference  has  a  purely  imaginary  spectrum  while  the  second  term  in  (6.3),  which  is 
just  a  (2r  —  1)*^  order  upwind  biased  difference,  has  a  spectrum  of  the  form 

2^^-^ Tr  ^sin  0  (sin  ^  +  i  cos  ^ )  (6.4) 

where  0  <  ^  <  2ir.  (6.3)  and  (6.4)  together  imply  that  the  spectrum  of  the  operator  L  lies 
fully  on  the  left  half  of  the  complex  plane.  Therefore,  with  an  appropriately  chosen  CFL 
number,  the  first  variation  of  the  WENO  schemes  are  /2-stable  when  the  3’’“^  or  higher  order 
Runge-Kutta  time  discretization  is  used. 

Let’s  define  by  u(xo,to,Ax)  the  numerical  solution  at  (xo,to)  €  R‘^  x  R+  for  grid  size 
Ax  and  fixed  CFL  number.  For  general  scalar  conservation  laws,  the  same  analysis  gives 

Theorem  6.1  For  the  initial  value  problem  of  (1.1)  with  m  =  1  (i.e.,  scalar  conservation 
laws),  V(xo,to)  €  R'^xR'*',  if  the  exact  solution  v  and  ^,g  have  r-|-[^^]-j-go+2  continuous 
derivatives  in  the  domain  of  dependence  of  (xo,to)  os  defined  in  [18],  the  WENO  schemes 
using  a  smooth  flux  splitting  and  a  order  Runge-Kutta  scheme  (n>  max(r,  3)  )  satisfy 

u(xo,  to,  Ax)  =  t;(xo,  to)  -I-  0{Ax'')  (6.5) 

for  appropriately  chosen  CFL  number.  Here  go  is  a  small  constant  integer  (see  [18]). 

For  a  few  special  cases,  we  list  the  CFL  numbers  in  Table  5. 

Table  5:  CFL  numbers  (n:  order  of  the  Runge-Kutta  scheme). 


n  =  3 

n  =  4 

r  =  2 

1.625 

1.745 

r  =  3 

1.434 

1.731 

7  Efficiency  Comparison 

In  this  section,  we  compare  the  efficiency  of  WENO-LF-4,  WENO-LF-5,  WENO-LF-5-PS 
and  ENO-LF-4  on  a  vector  supercomputer  (CRAY  C-90)  and  two  serial  workstations  (SUN 
SparclO  and  SGI  Indigo2). 
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ID,  2D  and  3D  Euler  systems  are  solved.  The  3D  Euler  system  is  (1.1)  with  d  =  3,  m  =  5 

and 

u  =  {p,pu,pv,pw,Ef-, 
f(u)  =  {pu,P  +  pu^,puv,puw,u(E  + 
g(u)  =  {pv,pvu,P  +  pv^,pvw,v(E  +  P))'^; 
h(u)  =  (p'w,pwu,pwv,P  +  pvP,w(E  +  P))^. 

where 

■P  =  (t  -  i){E  -  ^p(u^  +  v^  +  w^)) 

The  initial  condition  is 

/9  =  1  +  0.2  sin(7r(x  +  y  +  z)),  u  =  v  =  w  =  1,  P  =  1 

Here  we  use  f,g,h,x,y,z  instead  of  Xi,X2>®3-  The  ID  and  2D  Euler  systems  and 

their  initial  conditions  can  be  deduced  from  the  above  3D  problem  by  removing  the  extra 
degree  of  freedom(s). 

We  display  the  CPU  time  of  ENO-LF-4,  WENO-LF-4,  WENO-LF-5  and  WENO-LF- 
5-PS  (all  with  RK-4)  on  the  CRAY  C-90,  the  SparclO  and  the  SGI  Indigo2  in  Table  6. 
We  observe  that  WENO-LF-4  and  WENO-LF-5  are  at  least  twice  as  fast  as  ENO-LF-4  on 
the  CRAY  C-90  and  WENO-LF-5-PS  is  2.5  times  as  fast  as  ENO-LF-4  for  the  ID  Euler 
problem,  3.2  times  as  fast  for  the  2D  Euler  problem  and  3.9  times  as  fast  for  the  3D  Euler 
system.  On  the  workstations,  WENO-LF-4  and  WENO-LF-5  are  a  bit  faster  than  ENO- 
LF-4  on  the  SUN  SparclO  but  a  bit  slower  on  the  SGI  Indigo2.  WENO-LF-5-PS  is  1.5 
to  2.2  times  as  fast  as  ENO-LF-4  on  the  SUN  SparclO  and  on  the  SGI  Indigo2.  As  a 
reference,  we  also  include  the  CPU  times  of  a  typical  second  order  TVD  scheme  [8]  (Van 
Leer’s  limiter  with  2”"^  order  Runge-Kutta  scheme  in  time,  our  own  implementation)  in  the 
following  tables.  We  can  see  the  2^  order  scheme  is  about  10  times  as  fast  as  ENO-LF-4 
on  the  CRAY  C-90,  4.5  times  as  fast  on  the  SUN  SparclO  and  3.5  times  as  fast  on  the  SGI 
Indigo2. 

In  Table  7,  the  number  of  floating  point  operations  and  the  MFlops  (million  floating¬ 
point  operations  per  second)  are  given  for  the  2”*^  order  scheme,  ENO-Roe-4,  ENO-LF-4, 
WENO-LF-5,  WENO-LF-5-PS,  ENO-Roe-4-A  and  WENO-LF-5-A.  The  operation  count 
and  MFlops  for  WENO-LF-4  is  about  the  same  as  those  for  WENO-LF-5,  thus  omitted 
in  the  table.  We  can  see  all  the  WENO  schemes  achieve  the  speed  of  about  500  MFlops, 
which  is  50%  of  the  peak  speed  of  CRAY  C-90.  The  decrease  of  MFlops  for  high  dimensions 
is  because  of  the  shorter  array  length  N  used  in  our  tests.  Notice  also  that  the  operation 
count  per  grid  point  per  Runge-Kutta  stage,  of  the  full  characteristic  based  4*^^  or  5*^  order 
ENO  schemes  using  Lax- Friedrichs  building  blocks,  is  about  3  to  4  times  that  of  the  2’^'^ 
order  schemes.  This  ratio  actually  decreases  to  only  about  1.5  if  the  Roe  building  block 
is  used  instead,  i.e.  and  f~(u)  are  not  approximated  separately.  This  is  somewhat 

surprising,  as  it  was  commonly  believed  that  high  order  methods  are  much  more  expensive 
than  lower  order  ones.  When  the  Roe  building  block  is  used,  Yang’s  artificial  compression 
causes  a  40%  increase  in  operation  count  for  the  ID  Euler  system  and  a  65%  increase  for 
the  2D  Euler  system  as  we  can  see  from  the  operation  coimts  for  ENO-Roe-4  and  ENO- 
Roe-4-A.  When  the  Lax- Friedrichs  building  block  is  used,  the  increase  of  operation  count  is 
65%  for  the  ID  Euler  system  and  100%  for  the  2D  Euler  system  as  shown  by  the  operation 


(7.1) 

(7.2) 

(7.3) 

(7.4) 
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Table  6:  CPU  time  in  seconds.  N  points  in  each  spatial  dimension;  10"*“^  iterations  for  the 
d-dimensional  system. 


Q 

N 

ENO-LF-4 

WENO-LF-4 

WENO-LF-5 

WENO-LF-5-PS 

CRAY  C-90,  compiled  with  “-0  vectorS” 

D 

1.75 

16.67 

7.44 

7.45 

6.29 

B 

13.13 

122.52 

63.93 

60.84 

37.67 

60 

15.48 

171.42 

76.79 

78.89 

43.47 

SUN  SparclO  (66MHzj  HyperSparc),  compiled  with  -04” 

D 

1600 

69.43 

311.22 

317.55 

319.02 

215.95 

0 

2582.25 

2132.50 

2116.53 

1163.72 

40 

807.75 

716.05 

754.88 

389.77 

SGI  Indigo2  (75MHz,  R8000),  compiled  with  “-r8  -OS” 

D 

21.03 

66.21 

73.88 

77.01 

58.14 

0 

200 

151.26 

555.51 

564.54 

578.22 

347.48 

0 

167.44 

626.92 

699.58 

715.91 

366.29 

counts  for  WENO-LF-5  and  WENO-LF-5-A.  The  increase  in  CPU  time  is  well  reflected  by 
the  above  percentages. 

8  Numerical  Results 

8.1  Scalar  Conservation  Laws  in  One  Dimension 


Example  1.  Linear  Equation.  We  solve  the  linear  equation: 

Ut  +  Ux  =  0  -  1  <  X  <  1, 

u(x,0)  =  uo(x)  periodic 


where 


Uo(®) 


'  g  {Gix,  f3,z-S)  +  G(x,  l3,z  +  6)  +  4G(x,  I3,  z)) 
1 

^  l-|10(x- 0.1)1 

I  {F(x,  a,a-6)  +  F{x,  a,  a  +  ^)  +  4F(x,  a,  a)) 
(  0 


-0.8  <  X  <  -0.6; 

-0.4  <  X  <  -0.2; 
0  <  X  <  0.2; 

0.4  <  X  <  0.6; 
otherwise. 


G(x,^,z)  =  e-^(®-^)' 

F(x,  a,  a)  =  ^niax(l  —  a\x  -  o)2,0) 

The  constants  are  taken  as  c  =  0.5,  z  =  -0.7, 6  =  0.005,  a  =  10  and  /?  =  The  solution 
contains  a  smooth  but  narrow  combination  of  Gaussians,  a  square  wave,  a  sharp  triangle 
wave,  and  a  half  ellipse. 

We  compute  the  solution  up  to  t  =  8  with  200  points.  The  results  are  shown  in  Figure  2. 
We  observe  that  both  WENO- Roe-4  and  WENO-Roe-5  perform  better  than  ENO-Roe-4 
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Table  7:  Number  of  operations  per  Runge-Kutta  stage  per  grid  point  and  MFlops  on  CRAY 
C-90.  d:  the  spatial  dimension. 


Scheme 

Q 

100 

xly 

inai 

sign{x) 

X^,  y/x 

MFlops 

2"*^  order 

D 

82 

83 

9 

8 

3 

3 

478 

0 

239 

mm 

20 

8 

6 

400 

0 

476 

506 

39 

36 

15 

9 

mmnm 

ENO-Roe-4 

O 

WEm 

98 

3 

19 

0 

3 

0B0 

0 

309 

304 

6 

50 

0 

6 

191 

0 

663 

656 

9 

93 

0 

9 

— 

ENO-LF-4 

0 

244 

233 

3 

39 

0 

3 

1^00 

0 

6 

102 

0 

6 

IB0 

0 

1751 

1718 

9 

0 

9 

1001 

WENO-LF-5 

D 

235 

284 

27 

mm 

0 

3 

0 

703 

838 

70 

6 

0 

6 

503 

0 

129 

9 

0 

9 

I0E^H 

WENO-LF-5-PS 

0 

145 

129 

13 

3 

0 

4 

474 

0 

341 

315 

26 

6 

0 

8 

453 

0 

576 

39 

9 

0 

12 

I0S0 

ENO-Roe-4-A 

0 

144 

■EB 

4 

33 

6 

3 

0 

511 

484 

10 

106 

24 

6 

1001 

WENO-LF-5-A 

0 

375 

447 

37 

19 

12 

3 

0 

1379 

1654 

110 

70 

48 

6 

482 
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for  all  the  four  types  of  waves  in  the  initial  condition,  WEN 0- Roe-4  does  better  than 
WENO- Roe-5  at  acute  turns  in  the  solution  curve  (or  spike-like  peaks)  but  WENO-Roe-5 
does  better  for  the  square  wave  and  at  obtuse  turns  in  the  solution  curve.  With  Yang’s 
artificial  compression  technique,  WENO-Roe-5-A  performs  the  best  at  all  waves.  Note,  we 
have  adjusted  the  CFL  number  for  WENO-Roe-5-A  from  0.4  to  0.2.  For  CFL=0.4,  using 
a  =  40  in  (5.4)  gives  similar  results. 

Example  2.  N on-convex  Problems.  We  test  the  3’"'^  and  5*^  order  WENO  schemes  on 
the  Buckley- Leverett  problem  whose  flux  is 

^  4u2  +  (1  _  tt)2  (^-1) 

with  initial  data  w  =  1  in  [|,0]  and  u  =  0  elsewhere.  (For  the  numerical  results  of  the  4^^ 
order  WENO  scheme  of  Liu  et  al.,  see  [9]).  The  exact  solution  is  a  shock-rarefaction-contact 
discontinuity  mixture. 

The  results  obtained  by  WENO-RF-3  (with  RK-3)  and  WENO-RF^5  (with  RK-4)  are 
shown  in  Figure  3.  We  can  see  both  schemes  converge  to  the  correct  entropy  solution  and 
give  sharp  shock  profile.  Note  that,  around  discontinuities,  WENO  schemes  are  simulating 
the  base  ENO  schemes.  Therefore  the  sharpness  of  the  shock  profile  obtained  by  the  WENO 
schemes  are  only  expected  to  be  as  good  as  that  obtained  by  the  base  ENO  schemes. 
However,  in  terms  of  this  sharpness,  our  tests  show  that  the  3^*^  order  WENO  scheme  is 
comparable  to  the  3^^  order  ENO  scheme  instead  of  the  base  2^^  ENO  scheme  and,  the  5^^ 
WENO  scheme  is  comparable  to  the  4*^  order  ENO  scheme. 

8.2  Euler  System  in  One  Dimension 

Example  1.  ID  Riemann  Problems.  We  consider  here  two  well  known  problem,  which 
have  the  following  Riemann  type  initial  conditions: 

u(x,o)=(“'^ 

(  Ujt  if  X  >  0 

The  first  one  is  the  Sod  problem  [17].  The  initial  data  are: 

{PL,qL,PL)  =  (1,0,1);  {pR,qR,PR)  =  (0.125,0,0.1) 

The  second  one  is  the  Riemann  problem  proposed  by  Lax  [7]: 

{pL,  qi,  Pl)  =  (0.445, 0.698, 3.528);  (pR,  qR,  Pr)  =  (0.5, 0, 0.571) 

The  numerical  results  are  presented  in  Figure  4.  We  can  see  that  all  schemes  give  the 
correct  solution  with  good  resolution.  WENO-RF-5  is  better  than  WENO-LF-5-PS  which  is 
in  turn  better  than  WENO-RF-3.  We  note  that  Figures  4b  and  4c  (  Figures  4e  and  4f )  are 
comparable,  resp.  to  Figure  10  (  Figure  11 )  in  [15].  Also  see  Figures  9a  (  Figure  10a  )  ill  [9]- 
WENO-LF-5-A  does  much  better  than  all  other  schemes  at  the  contact  discontinuities.  We 
would  like  to  point  out  that,  according  to  our  experience  with  extensive  numerical  testing, 
these  two  problems,  especially  the  Lax’s  problem,  are  tough  test  cases  for  non- characteristic 
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Figure  2:  Linear  Equation.  Srd  order  Runge-Kutta  in  time.  200  points,  CFL=0.4  (0.2  for 
(d)  only),  T=8.  (a)  ENO-Roe-4.  (b)  WENO-Roe-4.  (c)  WENO-Roe-5.  (d)  WENO-Roe-5- 


t-0.4  cfl-0.8  n-ei  t-0.4  dl-0.8  n-81 


Figure  3:  The  Buckley-Leverett  problem,  (a);  WENO-RF-3;  (b):  WENO-RF-5. 

based  schemes  of  order  at  least  three.  Oscillations  can  easily  appear  for  such  schemes.  Here 
WENO-LF-5-PS  performs  well  in  these  two  cases. 

Example  2.  ID  Shock  Entropy  Wave  Interaction.  In  this  example,  we  test  the  WENO 
schemes  on  a  model  that  involves  a  moving  shock  interacting  with  an  entropy  wave  of  small 
amplitude.  On  a  domain  [0,5],  the  initial  condition  is: 

p  =  3.85714;  u  =  2.629369;  P  =  10.33333;  when  x  <  0.5 
«  =  0;  p  =  l;  when  a:  >  0.5 

where  e  and  k  are  the  amplitude  and  wave  number  of  the  entropy  wave,  resp.  The  mean 
flow  is  a  pure  right  moving  Mach  3  shock.  If  c  is  small  compared  to  the  shock  strength,  the 
shock  will  march  to  the  right  at  approximately  the  non-perturbed  shock  speed  and  generate 
a  sound  wave  which  travels  along  with  the  flow  behind  the  shock.  At  the  same  time,  the 
perturbing  entropy  wave,  after  “going  through”  the  shock,  is  compressed  and  amplified  and 
travels  approximately  at  the  speed  of  tt  +  c  where  u  and  c  are  the  velocity  and  speed  of  the 
sound  of  the  mean  flow  left  of  the  shock.  The  amplification  factor  for  the  entropy  wave  can 
be  obtained  by  linear  analysis.  See  [10,  21]  for  details.  In  order  to  get  rid  of  the  transient 
waves  due  to  the  non-numerical  initial  shock  profile,  we  let  the  shock  move  up  to  a;  =  4.5 
and  then  shuffle  it  back  to  a:  =  0.5.  The  solution  is  examined  when  the  shock  reaches 
X  =  4.5  the  second  time. 

The  goal  of  this  test  is  to  examine  the  stability  and  accuracy  of  the  WENO  schemes  in 
the  presence  of  the  shock.  Since  the  entropy  wave  here  is  set  to  be  very  weak  relative  to  the 
shock,  any  excessive  oscillation  could  pollute  the  generated  waves  (e.g.  the  sound  waves) 
and  the  amplified  entropy  waves.  In  our  tests,  we  take  e  =  0.01  and  k  —  13.  The  amplitude 
of  the  amplified  entropy  waves  predicted  by  the  linear  analysis  is  0.08690716  (shown  in 
the  following  figures  as  horizontal  solid  lines).  First  we  use  800  points  which  is  effectively 
20  points  in  each  wave  length  of  the  generated  entropy  wave.  Since  the  generated  sound 
waves  (or  pressure  wave)  are  of  lower  frequency  than  the  amplified  entropy  waves,  they  are 
much  better  resolved  by  this  grid  size.  Therefore  we  only  display  the  entropy  component 
of  the  numerical  solutions.  WENO-LF-4,  WENO-LF-5  and  WENO-LF-5-PS  are  used  in 
our  tests  and  the  results  are  shown  in  Figure  5  (the  mean  flow  has  been  subtracted  from 
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the  mimerical  solution).  We  see  that  all  three  schemes  catch  the  amplified  entropy  waves 
quite  well.  WENO-LF-4  performs  the  best  on  this  grid  and  WENO-LF-5  ranks  the  second. 
In  order  to  examine  the  relative  performance  of  WENO-LF-4  and  WENO-LF-5 ,  we  run 
the  same  test  on  a  grid  of  1200  points.  The  results  for  these  two  schemes  are  displayed  in 
Figure  6.  We  can  see  that  on  this  grid  (approximately  30  points  per  wave  length),  WENO- 
LF-5  is  as  accurate  as  WENO-LF-4.  In  fact,  on  finer  grids,  WENO-LF-5  becomes  more 
accurate  than  WENO-LF-4.  This  is  in  good  agreement  with  our  accuracy  test  in  section  4. 
For  the  purpose  of  comparison  with  low  order  schemes,  we  also  include  the  entropy  computed 
by  a  typical  second  order  scheme  [8]  (half  Van  Leer’s  limiter,  half  Superbee  limiter  with 
2"^^  order  Runge-Kutta  scheme  in  time,  2000  points).  The  advantage  of  using  higher  order 
schemes  for  this  example  is  apparent. 

Example  3.  Two  Interacting  Blast  Waves.  We  consider  here  the  interaction  of  two 
blast  waves.  The  initial  data  are  the  following: 

{U£,  if  0  <  a:  <  0.1 
UM  if  0.1  <  a;  <  0.9 
u/j  if  0.9  <  a:  <  1 

where 


PL  =  PM  =  Pfi  =  1  UL  =  UM  =  UR  =  0  Pl  =  10®  Pm  =  10  ®  Pr  =  10® 

A  reflective  boundary  condition  is  applied  at  both  a;  =  0  and  x  =  1.  See  [19]  for  a  detailed 
discussion  of  this  problem. 

Three  grids  are  used:  199,  399,  799  points.  We  examine  our  numerical  solutions  at 
t  =  0.038.  The  “exact”  solution  (solid  lines  in  all  the  pictures)  are  computed  by  ENO-RF-5 
with  1600  points.  In  Figure  7,  we  show  the  density  computed  by  WENO-RF-3  (with  RK-3), 
WENO-RF-4,  WENO-RF-5  and  WENO-RF-5-A  (with  RK-4). 

We  observe  that  the  4‘^  order  and  order  WENO  schemes  are  much  better  than 
the  3”'^  order  WENO  scheme  and  the  results  are  comparable  with  those  obtained  by  the 
unmodified  ENO-RF-4  (see  Figure  12  in  [15].  Note,  the  4*^  order  ENO  scheme  in  the 
norm  was  denoted  as  ENO-RF-3  there).  WENO-RF-4  is  slightly  better  than  WENO-RF-5 
on  the  medium  grid  while  on  the  fine  grid  WENO-RF-5  seems  to  be  better.  The  results 
of  WENO-RF-5- A  on  the  coarse  and  medium  grids  are  nearly  as  good  as  WENO-RF-5  on, 
resp.  medium  and  fine  grids. 

Example  4-  Quasi-One  Dimensional  Nozzle  Flow.  In  this  example,  we  use  the  WENO 
schemes  to  solve  the  steady  state  quasi-lD  nozzle  flow.  The  governing  equation  of  the 
quasi-lD  nozzle  flow  is  the  ID  Euler  system  with  the  following  forcing  term: 

g(u,x)  =  -^{pu,pu'^,u{E-\- P)f 

where  A  =  A{x)  is  the  cross  area  function  of  the  nozzle  and  Ax  =  ^.  The  nozzle  here  is 
of  unit  length,  whose  shape  is  determined  by  assuming  a  linear,  isentropic  Mach  number 
distribution,  which  is  0.8  at  x  =  0  (the  entrance)  and  1.8  at  x  =  1  (the  exit).  The  exit  flow 
condition  is  then  decided  by  the  prescribed  shock  position,  which  is  x  =  0.5  in  our  test. 

In  Figure  8,  we  display  the  density  computed  by  WENO-Roe-4  and  WENO-Roe-5  with 
34  points.  We  can  see  both  schemes  converge  nicely  to  the  exact  solution  (solid  line  in  the 
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ENTROPY  WENO-LF-5/RK-4  N=800 


ENTROPY  WENO-LF-5-PS/RK-4  N=800 


0  1 

2  3 

ENTROPY 

LAX-LIU/RK-2  N=2000 

Figure  5:  ID  Shock  entropy  wave  interaction.  Entropy.  At  =  0.6Ax 


ENTROPY  WENO-LF-4/RK-4  N=1200 


ENTROPY  WENO-LF-5/RK-4  N=1200 


Figure  6:  ID  Shock  entropy  wave  interaction(cont’d).  Entropy.  At  =  0.6Aa:. 

pictures).  The  residue  computed  with  both  schemes  settles  down  to  10“^  for  this  grid,  and 
to  a  smaller  number  for  a  finer  grid. 

This  example  shows  that  WENO  has  its  advantage  for  steady  state  calculations. 

8.3  Euler  System  in  Two  Dimensions 


Example  1.  Oblique  Sod’s  Problem.  The  purpose  of  this  test  is  to  analyze  the  capability 
of  WENO  schemes  in  resolving  waves  that  are  oblique  to  the  computational  mesh.  The 
2D  Sod’s  problem  is  solved  where  the  initial  jump  makes  an  angle  6  against  the  a:-axis 
(0  <  ^  <  f ).  K  ^  we  have  the  one-dimensional  Sod’s  problem.  If  0  <  0  <  f ,  aU 
the  waves  produced  will  be  oblique  to  the  rectangular  computational  mesh.  We  take  our 
computational  domain  to  be  [0,6]  x  [0, 1]  and  position  the  initial  jump  at  (®,y)  =  (2.25,0). 
The  physical  domain  varies  with  0  and  is  taken  as  [0,;j|g]  x  [0,51^].  The  scaling  factor 
is  to  ensure  the  same  grid  resolution  normal  to  the  wave  propagation  on  a  given  mesh 
at  some  fixed  time  for  all  choices  of  6.  See  [3]  for  details.  We  take  0  to  be  arctan  1,  arctan2, 
and  arctan  4.  The  solution  is  computed  up  to  t  =  1.2  on  a  96  x  16  mesh. 

WENO-LF-4,  WENO-LF-5,  WENO-LF-5-A  and  WENO-LF-5-PS  are  used  in  our  tests 
and  the  results  are  compared  with  that  obtained  by  ENO-LF-4  (all  with  RK-4  and  At  = 
O.GAx).  For  the  case  0  =  arctan  1,  we  display  the  density  contours  obtained  by  WENO- 


25 


DENSITY 


t»0.03a  cfl.0.6  n-399 


DENSITY 


t.0.038  cf]»0.6  n.799 


Figure  7:  Two  interacting  blast  waves 


Figure  8:  Density.  Steady  quasi-lD  nozzle  flow.  34  points.  RK-3  in  time. 

LF-5-PS  in  Figure  9a;  In  Figure  9b,  we  show  the  densities  at  y  =  0  obtained  by  all  four 
schemes.  We  can  see  that  all  WENO  schemes  are  doing  well  in  resolving  the  oblique  waves 
and  their  differences  from  the  ENO-LF-4  (except  WENO-LF-5-A)  are  barely  noticeable. 
WEN0-LF-5-A  gives  sharp  profile  of  the  contact  discontinuity  as  expected.  In  Figure  9(c- 
f),  a  more  quantitative  study  is  carried  out.  Namely,  for  each  scheme,  we  measure  the 
differences  between  oblique  cases  and  the  one-dimensional  case.  We  can  see  that  WENO- 
LF-4  and  WENO-LF-5  perform  similarly  as  ENO-LF-4  does  whUe  WENO-LF-5-PS  gives 
a  slightly  larger  deviation  near  the  contact  discontinuity.  However,  this  small  difference 
can  be  regarded  as  negligible.  We  want  to  note  that  WENO-LF-5-PS  performs  weU  at  the 
shock. 

Example  2.  A  Mach  3  Wind  Tunnel  with  a  Step.  This  model  problem  has  been  carefully 
examined  in  [19].  The  set  up  of  the  problem  is  the  following:  The  wind  tunnel  is  1  length 
unit  wide  and  3  length  units  long.  The  step  is  0.2  length  units  high  and  is  located  0.6  length 
units  from  the  left-hand  end  of  the  tunnel.  The  problem  is  initialized  by  a  right-going  Mach 
3  flow.  Reflective  boundary  conditions  are  applied  along  the  walls  of  the  tunnel  and  an  in¬ 
flow  and  an  out-flow  boundary  conditions  are  applied  at  the  entrance  (left-hand  end)  and 
the  exit  (right-hand  end).  For  the  treatment  of  the  singularity  at  the  corner  of  the  step,  we 
adopt  the  same  technique  used  in  [19],  which  is  based  on  the  assumption  of  a  nearly  steady 
flow  in  the  region  near  the  comer. 

WENO-LF-4  and  WENO-LF-5  are  used  in  our  tests  and  the  results  are  compared  with 
those  obtained  by  ENO-LF-4  (all  with  RK-4  and  At  =  0.6Ax).  Two  grids  are  used:  122  x  39 
and  242  x  79.  They  correspond  resp.  to  the  medium  and  fine  grids  in  [19]. 

In  Figure  10  to  11,  we  show  the  density  component  obtained  by  all  three  schemes  on  the 
two  grids.  We  can  see  that  aU  the  schemes  perform  well  with  good  resolution.  Relatively 
speaking,  WENO-LF-4  and  WENO-LF-5  have  slightly  better  resolution  at  the  contact  line 
(originated  from  the  Mach  step)  and  contain  less  visible  “bumps”,  which  are  indeed  small 
numerical  oscillations,  than  ENO-LF-4. 

WENO-LF-5-PS  does  not  work  for  this  problem  because  of  the  strong  reflecting  waves. 
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Figure  9:  Oblique  Sod’s  problem,  (a)  Density  contours.  WENO-LF-5-PS,  9  =  arctanl.  (b) 
Density,  y  =  0,  ^  =  arctan  1.  (c)-(f)  pe- Rid,  y  =  0.  (c)  ENO-LF-4.  (d)  WENO-LF-4.  (e) 
WENO-LF-5.  (f)  WENO-LF-5-PS. 
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Figure  11:  Flow  past  a  forward  facing  step  (conf 
WENO-LF-4.  (b)  WENO-LF-5.  (c)  ENO-LF-4. 


cfl=0.6  Grid:  242x79  t=4. 


WENO-LF-5 


ENO-LF-4 


Example  3.  Double  Mach  Reflection  of  a  Strong  Shock.  The  computational  domain 
for  this  problem  is  chosen  to  be  [0,4]  x  [0,1].  The  reflecting  wall  lies  at  the  bottom  of 
the  computational  domain  starting  from  x  —  Initially  a  right-moving  Mach  10  shock  is 
positioned  at  x  =  |,  j/  =  0  and  makes  a  60°  angle  with  the  x-axis.  For  the  bottom  boundary, 
the  exact  post-shock  condition  is  imposed  for  the  part  from  x  =  0  to  x  =  |  and  a  reflective 
boundary  condition  is  used  for  the  rest.  At  the  top  boundary  of  our  computational  domain, 
the  flow  values  are  set  to  describe  the  exact  motion  of  the  Mach  10  shock.  See  [19]  for  a 
detailed  description  of  this  problem. 

Two  grids  have  been  used  in  our  tests:  240  x  59  and  480  x  119.  They  correspond  to  the 
medium  and  fine  grids  in  [19],  resp.  We  will  only  show  the  solutions  on  part  of  the  domain: 
[0, 3]  X  [0, 1]  where  most  of  the  flow  features  are  located. 

We  use  WENO-LF-4,  WENO-LF-5  and  ENO-LF-4  (all  with  RK-3  and  At  =  0.6Ax)  in 
our  tests.  We  show  the  density  contours  obtained  by  these  three  schemes.  See  Figure  12 
and  13.  We  see  that  aU  three  schemes  resolve  the  two  Mach  stems  well.  Again  WENO-LF- 
5-PS  does  not  work  because  of  the  strong  reflecting  wave  pattern  in  this  problem. 

Example  4-  2D  Shock  Entropy  Wave  Interaction.  In  this  example,  we  test  the  WENO 
schemes  on  a  2D  model  that  involves  the  interaction  between  a  normal  shock  and  a  weak 
entropy  wave  which  makes  an  angle  9t  €  (0,7r/2)  against  the  x-axis.  If  =  0,  we  have 
essentially  the  ID  problem  (see  Example  2  in  Section  8.2).  Since  the  weak  entropy  waves 
are  oblique  to  the  shock,  the  waves  generated  by  the  interaction  are  much  more  difficult 
to  resolve  than  in  the  ID  case.  Our  goal  here  is  to  further  examine  the  capability  of  the 
WENO  scheme  in  capturing  such  small  scale  waves  in  the  presence  of  a  shock.  See  [21,  16] 
for  detailed  discussions  on  this  subject.  The  set-up  of  the  problem  is  the  following:  for  a 
right  moving  normal  shock  of  Mach  number  M,  we  add  a  small  entropy  wave  to  the  flow 
on  the  right  of  the  shock  which  is  equivalent  to  changing  only  the  density  of  the  flow  on  the 
right  of  the  shock  to: 

p  =  p^e~^r(smpr)IPr 

where  pr  and  Pr  are  resp.  the  density  and  pressure  of  the  right  state  of  the  shock,  = 
kr{x  cos  Ot  -H  y  sin  9r)  and  kr  is  the  entropy  wave  number.  In  order  to  enforce  periodic 
boundary  conditions  in  the  j/-direction,  we  take  the  computation  domain  to  be  [0,5]  x 
[0,  XTsme^]-  initially  position  the  normal  shock  at  x  =  0.5  and  allow  it  to  move  up  to 
X  =  4.5  and  then  shuffie  the  data  back  to  x  =  0.5.  We  extract  the  data  at  the  time  when 
the  shock  moves  up  to  x  =  4.5  again.  See  [16]  for  a  similar  implementation. 

In  our  tests,  we  take  Ur,  the  velocity  on  the  right  of  the  shock,  to  be  0  and  set  M  =  Z, 
€r  =  0.01,  fcr  =  15,0,.  =  30°. 

We  measure  the  performance  by  comparing  the  amplitude  of  the  amplified  entropy 
waves,  which  is  computed  by  a  Fourier  analysis  in  the  y  direction  for  all  fixed  grid  values 
X  €  [3.4, 4.4]. 

Both  ENO-LF-4  and  ENO-LF-5  suffer  a  loss  of  accuracy  if  not  modified  (see  [16]). 
WENO-LF-4  and  WENO-LF-5  work  nicely  without  any  modification. 

The  loss  of  accuracy  of  the  ENO  schemes  can  be  easily  fixed  by  the  techniques  introduced 
in  [1,  13],  which  effectively  force  an  upstream  centered  stencil  to  be  chosen  away  from  the 
shock  and  free  adapted  stencil  to  be  used  near  the  shock.  The  techniques  can  also  be  adapted 
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(a) 

1.0 

0.5 

0 

(b) 

1,0 

0.5 

0 

(c) 

Figure  13:  Double  Mach  reflection  (cont’d).  Density  on  fine  grid:  480  x  119.  (a)  WENO- 
LF-4.  (b)  WENO-LF-5.  (c)  ENO-LF-4. 


DENSITY  ENO-LF-4 


30  contours  from  1.731  to  20.92  Grid:  480x1 19  cfl=0.6  t=0.2 


DENSITY  WENO-LF-5 


30  contours  from  1.731  to  20.92  Grid:  480x1 19  cfl=0.6  t=0.2 


30  contours  from  1.731  to  20.92  Grid:  480x1 19  cfl=0.6  t=0.2 
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to  enhance  the  performance  of  WENO  schemes  by  modifying  the  weights  as  follows: 


Qk  = 


Cl  if  ISi  <  a  ^  for  any  /  =  0, . . . ,  r- 1, 
ojk  otherwise. 


(8.2) 


where  a  is  taken  as  2  and  =  0,...,r-l  are  the  regularly  computed  weights.  (8.2) 
leads  to  the  optimal  weights  being  used  for  stencils  away  from  the  shock  and  regularly 
computed  weights  being  used  near  the  shock.  We  denote  the  modified  ENO-LF-4,  ENO- 
LF-5,  WENO-LF-4  and  WENO-LF-5  to  be,  resp.  EN02-LF-4,  EN02-LF-5,  WEN02-LF-4 
and  WEN02-LF-5.  Note  that,  with  the  modified  weights,  WENO-LF-4  becomes  5*^  order 
accurate  in  smooth  regions.  In  all  tests,  RK-4  is  used  with  At  =  O.fiAx. 

First  we  use  800  points  in  the  x- direction  and  20  points  in  the  y-direction  which 
give  approximately  20  points  per  entropy  wave  length  in  both  directions.  In  Figure  14, 
the  amplitude  of  the  amplified  entropy  waves  obtained  by  all  aforementioned  schemes  are 
displayed.  The  solid  horizontal  line  is  the  amplitude  predicted  by  linear  analysis  which  is 
0.08744786  (see  [10,  21]).  We  see  that  the  modified  schemes  generally  perform  better  in 
terms  of  accuracy  and  decay  rate  than  the  unmodified  schemes.  As  a  reference,  we  have  also 
included  the  amplitudes  obtained  by  a  typical  second  order  TVD  scheme  [8]  (half  Van  Leer’s 
hmiter,  half  Superbee  limiter  with  2”^  order  Runge-Kutta  scheme  in  time,  800  points).  We 
can  see  that  high  order  schemes  perform  much  better  than  the  second  order  schemes  in 
terms  of  accuracy  and  decay  rate  for  this  problem. 

WENO-LF-5-PS  does  not  perform  well  for  this  problem  even  with  the  remedy  above. 
This  indicates  that  the  pressure-entropy  combination  is  not  good  enough  to  indicate  pre¬ 
cisely  the  smoothness  of  the  numerical  solution.  This  causes  oscillations  generated  at  shocks 
and  thus  destroys  the  accuracy  of  the  scheme  in  resolving  the  waves  which  have  “undergone” 
the  interaction  with  the  shock. 


Remark:  We  have  seen  that  WENO-LF-5-PS  does  not  work  for  the  step  problem 
and  the  double  Mach  reflection  problem  because  it  can  not  handle  the  reflective  boundary 
properly.  This  can  be  explained  by  the  following:  the  usual  way  of  imposing  the  reflective 
boundary  condition^  is  to  reverse  the  normal  velocity  at  the  grid  points  which  are  symmetric 
with  respect  to  that  boundary  while  setting  other  flow  quantities  (density,  pressure  and 
tangential  velocity)  to  be  the  same;  in  particular  the  pressure  and  entropy  at  each  pair  of 
symmetric  grid  points  are  identical.  Therefore  neither  the  pressure  nor  the  entropy  can 
indicate  possible  jumps  in  the  normal  velocity.  This  failure  will  result  in  an  unstable  weight 
distribution  in  the  normal  direction  near  the  reflective  boundary  and  cause  fatal  errors  such 
as  density  becoming  negative.  An  immediate  “fix”  seems  to  be  using  the  normal  velocity 
to  compute  the  weights  for  one  of  the  linearly  degenerate  fields.  Unfortunately,  the  jump 
in  the  normal  velocity  is  not  like  a  contact  discontinuity,  which  belongs  solely  to  one  of  the 
characteristic  fields.  While  this  might  cure  the  ill  distribution  of  the  weights  in  the  field,  in 
which  the  velocity  is  used,  it  can  not  cure  this  iU  distribution  in  other  fields. 

However,  WENO-LF-5-PS  can  be  applied  to  problems  where  the  reflective  boundary 
does  not  play  a  vital  role.  As  an  example,  we  look  at  the  following  model  problem. 

Example  5.  Shock  Vortex  Interaction.  This  model  problem  describes  the  interaction 
between  a  stationary  shock  and  a  vortex.  The  computational  domain  is  taken  to  be  [0, 2]  x 

*We  assume  here  the  physical  boundary  is  flat,  as  is  the  case  in  aforementioned  problems 
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Figure  14:  2D  shock  entropy  wave  interaction.  Amplitude  of  amplified  entropy  waves.  800 
points  (about  20  points  per  entropy  wave  length). 

[0, 1].  A  stationary  Mach  1.1  shock  is  positioned  at  a:  =  0.5  and  normal  to  the  x-axis.  Its 
left  state  is  {p,u,v,P)  =  (1,.^,0, 1).  A  small  vortex  is  superposed  to  the  flow  left  to  the 
shock  and  centers  at  {xc,yc)  =  (0.25,0.5).  We  describe  the  vortex  as  a  perturbation  to  the 
velocity  (n,  v),  temperature  (T  =  y )  and  entropy  (5  =  In  ^)  of  the  mean  flow  and  denote 
it  by  the  tilde  values: 


u  — 

sin  0 

(8.3) 

V  = 

— cos  9 

(8.4) 

f  = 

(7  — 

407 

(8.5) 

II 

0 

(8.6) 

where  t  =  ^  and  r  =  \/(x  -  «c)^  +  (y  —  t/c)^-  Here  e  indicates  the  strength  of  the  vortex, 
a  controls  the  decay  rate  of  the  vortex  and  Tc  is  the  critical  radius  for  which  the  vortex 
has  the  maximum  strength.  In  our  tests,  we  choose  c  =  0.3,  Tc  =  0.05  and  a  =  0.204.  The 
above  defined  vortex  is  a  steady  state  solution  to  the  2D  Euler  equation. 

We  use  a  grid  of  251  x  100  which  is  uniform  in  y  but  refined  in  x  around  the  shock  using 
a  Roberts  transform  (see  [2]  and  the  references  there).  The  upper  and  lower  boundaries 
are  intentionally  set  to  be  reflective.  The  pressure  contours  obtained  by  WENO-LF-5-PS 
at  t  =  0.05,  t  =  0.20  and  t  =  0.35  are  shown  in  Figure  15(a-c).  We  can  see  that  for  this 
problem,  where  the  reflective  boundary  is  nonessential,  WENO-LF-5-PS  works  nicely.  To 
appreciate  this  further,  we  look  at  the  solution  at  t  =  0.8.  By  this  time  one  branch  of  the 
shock  bifurcations  has  reached  the  top  boundary  and  been  reflected.  The  pressure  contours 
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obtained  by  WENO-LF-5-PS  at  this  moment  are  shown  in  Figure  15d.  We  see  that  the 
reflection  is  well  captured. 

In  Figure  15(e-g),  we  compare  the  results  obtained  by  WENO-LF-5-PS,  WENO-LF- 
5  and  ENO-LF-4.  90  contours  are  drawn  for  the  pressure  component  in  the  range  of 
(1.19,1.37).  We  see  that  the  three  methods  give  approximately  the  same  resolution.  A 
careful  examination  reveals  that  WENO  schemes  are  slightly  better  in  the  sense  that  less 
numerical  noise  is  generated.  Between  the  two  WENO  schemes,  WENO-LF-5  seems  a  little 
better  for  the  same  reason.  For  a  qualitative  comparison,  see  also  [2].  Note  that  a  different 
vortex  is  used  there. 


Figure  15:  2D  shock  vortex  interaction.  Pressure,  (a)-(d)  WENO-LF-5-PS.  30  contours, 
(a)  t=0.05.  (b)  t=0.20.  (c)  t=0.35.  (d)  t=0.80.  (e)-(g)  t=0.60.  90  contours  from  1.19  to 
1.37.  (e)  WENO-LF-5-PS.  (f)  WENO-LF-5.  (g)  ENO-LF-4. 


Example  6.  Flow  Past  a  Cylinder.  In  this  test,  we  use  the  WENO  schemes  to  simiilate 
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the  supersonic  flow  past  a  cylinder.  In  the  physical  space,  a  cylinder  of  unit  radius  is 
positioned  at  the  origin  on  a  x  —  y  plane.  The  computational  domain  is  chosen  to  be 
[0,1]  X  [0,1]  on  the  plane.  The  mapping  between  the  computational  domain  and  the 
physical  domain  is: 

X  =  (Rx  -  {Rx  -  IX)  cos{0{2i]  -  1))  (8.7) 

y  =  {Ry-{Ry-  1)0  sin((9(277  -  1))  (8.8) 

where  we  take  Rx  =  3,Ry  =  6  and  ^  •  See  [16]  for  the  eigenvalues  and  eigenvectors  of  2D 

Euler  systems  on  general  structured  grids.  A  uniform  mesh  of  60  x  80  in  the  computational 
domain  is  used.  For  an  illustration  of  the  mesh  in  the  physical  space  (drawing  every  other 
grid  hne),  see  Figure  16a. 

The  problem  is  initialized  by  a  Mach  3  shock  moving  toward  the  cylinder  from  the 
left.  A  reflective  boundary  condition  is  imposed  at  the  surface  of  the  cylinder,  i.e.  ^  =  1, 
inflow  boundary  condition  is  applied  at  ^  =  0  and  outflow  boundary  condition  is  applied  at 
7?  =  0,1, 

The  pressure  contour  obtained  by  WENO-LF-5  with  RK-4  and  At  =  0.6 A»  is  shown  in 
Figure  16b.  Similar  results  can  be  obtained  by  WENO-LF-4  and  ENO-LF-4. 

9  Conclusion 

With  the  new  smooth  measurement,  which  is  based  on  minimizing  the  norm  of  the 
derivatives  of  the  interpolation  polynomials,  the  WENO  schemes  formulated  from  the 
order  ENO  schemes  can  be  made  (2r— 1)*^  order  accurate  in  smooth  regions  of  the  flux 
function  (in  spatial  variables),  at  least  for  r  =  2,3.  However,  at  discontinuities,  all  WENO 
schemes  are  just  accurate  (r  is  the  order  of  the  base  ENO  scheme). 

The  4*^  order  WENO  scheme  of  Liu  et  al.  and  the  5*^  order  WENO  scheme  resulting 
from  the  new  smoothness  measurement  are  found  to  be  at  least  twice  as  fast  as  the  4‘^  order 
ENO  schemes  on  vector  supercomputers  (e.g.  CRAY  C-90)  and  as  fast  on  serial  machines 
(therefore  on  parallel  machines  as  well).  Many  ID  and  2D  numerical  tests  suggest  that  both 
WENO  schemes  are  very  robust  for  shock  calculations.  The  4*^  order  WENO  scheme  of 
Liu  et  al.  is  slightly  more  accurate  than  the  5*^  order  WENO  scheme  on  coarse  grids  (20 
points  or  less  per  wave  length)  but  becomes  less  accurate  on  finer  grids. 

For  Euler  systems,  we  also  suggest  computing  the  weights  from  pressure  and  entropy 
instead  of  the  projected  values.  The  restdting  WENO  schemes  are  about  twice  as  fast  as 
the  WENO  schemes  which  use  the  projected  values  to  compute  weights,  and  work  weU  for 
problems  which  do  not  contain  strong  shocks  or  strong  reflected  waves. 

More  detailed  numerical  results  for  WENO  schemes  can  be  found  in  [6]. 

We  have  also  adopted  the  artificial  compression  method  of  Yang  [20]  to  enhance  the  per¬ 
formance  of  WENO  schemes  at  contact  discontinuities.  However,  the  CPU  cost  is  increased 
by  as  much  as  100%  when  a  Lax-Friedrichs  biiilding  block  is  used.  We  believe  the  idea  of 
artificial  compression  method  can  be  adapted  directly  into  the  weight  definition  to  achieve 
the  sharpening  effect  at  a  much  lower  expense.  This  will  be  investigated  in  the  future. 
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PHYSICAL  GRID  30x40 


PRESSURE  60x80 
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Figure  16:  Flow  past  a  cylinder,  (a)  Physical  grid,  (b)  Pressure.  WENO-LF-5  with  RK-4. 
20  contours 


suggestions. 
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